DISCONTINUOUS GALERKIN METHOD FOR FRACTIONAL 
CONVECTION-DIFFUSION EQUATIONS 

Q. XU* AND J. S. HESTHAVENt 

Abstract. We propose a discontinuous Galerkin method for convection-subdiffusion equations 
with a fractional operator of order «(1 < a < 2) defined through the fractional Laplacian. The 

£f*l , fractional operator of order a is expressed as a composite of first order derivatives and fractional 

integrals of order 2 — a, and the fractional convection-diffusion problem is expressed as a system of 
low order differential/integral equations and a local discontinuous Galerkin method scheme is derived 

£vj ■ for the equations. We prove stability and optimal order of convergence 0(h kJtl ) for subdiffusion, 

and an order of convergence of 0(h + 2 ) is established for the general fractional convection-diffusion 
O 1 . problem. The analysis is confirmed by numerical examples. 
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1. Introduction. In this paper we consider the following convection diffusion 
equation 

^ + lf(u) = e(-(-Ar/2)u(x,t), xeR,te(0,T], (11) 

u(x,0) = u (x), i£l 
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where / is assumed to be Lipschitz continuous and the subdiffusion is defined through 
the fractional Laplacian — (— A) Q / 2 (a e (0,2]), which can be defined using Fourier 
analysis as[22j GES EH El , 

o 

^; - (^Ap/%(0 = (27r)°|ertt(0 (1-2) 

© ; Equivalent^, the fractional Laplacian can also be defined by a singular integral^, 

-(-A)*«(*) = «./ U{x+ .^ x) dz (1.3) 

J|z|>0 \ z \ 



Equation (jl.ip , also known as a fractional conservation law, can be viewed as a gener- 
alization of the classical convection diffusion equation. During the last decade, it has 
arisen as a suitable model in many application areas, such as geomorphology[25, 26, 2 , 
overdriven detonations in gases [HI [I], signal processing [5], anomalous diffusion in 
semiconductor growth [3B] etc. With the special choice of f(u) = w 2 /2, it is recog- 
nized as a fractional version of the viscous Burgers' equation. Fractional conservation 
laws, especially Burgers' equation, has been studied by many authors from a theo- 
retical perspective, mainly involving questions of well-posedness and regularity of the 
solutions [2 71 131 I3"T1 [BUI |2"U] . In the case of a < 1, the solution is in general not smooth 
and shocks may appear even with smooth initial datum. Similar to the classical scalar 



'School of Mathematics and Statistics, Central South University,410083 Changsha, China;Division 
of Applied Mathematics, Brown University, RI, 02912 USA. Email: qw_xuahotmail.com 

^Division of Applied Mathematics, Brown University, RI, 02912 USA. Email: 
Jan . HesthavenOBrown . edu 



2 Q. Xu and J. S. Hesthavon 

conservation laws, an appropriate entropy formulation is needed to guarantee well- 
posedness[21H]- For the case a > 1, the existence and uniqueness of a regular solution 
has been proved in [3J [3U] • In this case, the nonlocal term serves as subdiffusion term, 
and smooth out discontinuities from initial datum. 

Numerical studies of partial differential equations with nonlocal operator have 
attracted a lot of interest in recent years. Liu, Deng et al. have worked on numerical 
methods for fractional diffusion problems with fractional Laplacian operators or Riesz 
fractional derivatives [17l |40j [29j [19] . Briani, Cont, Matache, et al. have considered 
numerical methods for financial model with fractional Laplacian operators P [TBI 152"] . 
However, for conservation laws with fractional Laplacian operators, the development 
of accurate and robust numerical methods remains limited. Droniou is the first to 
analyze a general class of difference methods for fractional conservation laws [21 . Az- 
eras proposed a class of finite difference schemes for solving a fractional antidiffusive 
equation[B]. Bouharguane proposed a splitting methods for the nonlocal Fowler equa- 
tion in [Sj- For the case a < 1, Cifani et al. [THEI] applied the discontinuous Galerkin 
method to the fractional conservation law and degenerate convection diffusion equa- 
tions and developed some error estimation. Unfortunately, their numerical results 
failed to confirm their analysis. 

The discontinuous Galerkin method is a well established method for classical con- 
servation laws |2H]. However, for equations containing higher order spatial derivatives, 
discontinuous Galerkin methods cannot be directly applied [15l 39 . A careless appli- 
cation of the discontinuous Galerkin method to a problem with high order derivatives 
could yield an inconsistent method. The idea of local discontinuous Galerkin meth- 
ods [15) for time dependent partial differential equations with higher derivatives is 
to rewrite the equation into a first order system and then apply the discontinuous 
Galerkin method to the system. A key ingredient for the success of this methods is 
the correct design of interface numerical fluxes. These fluxes must be designed to 
guarantee stability and local solvability of all the auxiliary variables introduced to 
approximate the derivatives of the solution. 

In this paper, we will consider fractional convection diffusion equations with a 
fractional Laplacian operator of order a(l < a < 2). Especially for 1 < a < 2, it is 
conceptionally similar to a fractional derivative with an order between 1 and 2. To 
obtain a consistent and high accuracy method for this problem, we shall rewrite the 
fractional operator as a composite of first order derivatives and a fractional integral, 
and convert the fractional convection diffusion equation into a system of low order 
equations. This allows us to apply the local discontinuous Galerkin method. 

This paper is organized as follows. In Section 2, we introduce some basic defi- 
nitions and state a few central lemmas. In Section 3, we will derive the discontin- 
uous Galerkin formulation for the fractional convection-diffusion law. In section 4, 
we develop stability and convergence analysis for fractional diffusion and convection- 
diffusion equations. In Section 5 some numerical examples are carried out to support 
our analysis. Conclusions are offered in Section 6. 

2. Background and definitions. Apart from the definitions of the fractional 
Laplacian based on the Fourier and integral form above, it can also be defined using 
ideas of fractional calculus [2"2l [551 |4T)] . as 
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where _ oc D" and ^-DJ, refer to the left and right Riemann-Liouville fractional deriva- 
tives, respectively, of a'th order. This definition is also known as a Riesz derivative. 
In this paper, we will based our developments and analysis mainly on this defini- 
tion and will, in preparation, introduce a few definitions and properties of fractional 
integral/derivative to set the stage. 

The left and right fractional integrals of order a are defined as 

-ooI>(x) = -jL T (* - t) a - l u{t)dt (2.2) 



x I^u(x) = = n; (t- x) a - l u(t)dt (2.3) 

r (") Jx 

This allows for the definition of the left and right Riemann-Liouville fractional deriva- 
tive of order a (n — 1 < a < n) as, 

_ 00 D«u(x) = - — / (x-t) n - l - a u(t)dt = D n (_ 00 I^u(x)) (2.4) 

T{n-a) \dxj J _ 00 

1 f d \ n f 00 
x D^u(x) = (--) J (t xY- l -*u{t)dt = {-DY{J^*u{x)) (2.5) 

Mirroring this, Caputo's left and right fractional derivatives of order a (n — 1 < 
a < n) arc defined as 

_gD2u{x) = *_ - J X Jx ~ tr- 1 -" (j^ u{t)dt = _^-(0\(x)) (2.6) 

C x D^u(x) = T{n _ a) J (t-xT- 1 - (--J u{t)dt = x IZ~ a ((-D) n u(x)) (2-7) 

The two definitions, the Riemann-Liouville fractional derivative and Caputo's frac- 
tional derivatives, are naturally related and they are equivalent provided all derivatives 
less than order n of u(x) disappear when x — > ±oo. 

Fractional integrals and derivatives satisfy the following properties, 

Lemma 2.1. (Linearity 1341) ■ 

_ oc i:(\f(x)+ f ig(x)) = \_ 00 r:f(x)+ f i_ 00 I°g(x) (2.8) 

^D^Xfix) + W (s)) - A _„,!£/(») + /i . OB D^g{x) (2.9) 

Lemma 2.2. (Semigroup property [34 1) ■ For a, (3 > 

^/^/(Z) = _ 00 4 Q (-oo-tf/O*)) = -oo-tf (-«,-£/(*)) (2T0) 

Lemma 2.3. /5^/ Suppose u^(x) — O.when x — >• ±oo, V0 ^ j ^ n ( n— I < a < 
n), then 

_ 0o D>(x) = £>" (^/r^W) = .^-" (!>"«(*)) (2.11) 

„££«(*) - (-£)" ( B ir a u(a ! )) = x /r Q ((--D) n «(a:)) (2.12) 
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From Lemma 12.31 we get, for 1 < a < 2 and a continuous function u with 
u V\x) = 0, £ -> ±oo, ^ j ^ 2 



1 d 



\- Xx dx +Xo ° dx 



2cos(a7r/2) dx 

1 /" r2 _ a d 2 n(a?) , 2 _ a rfVx) 

— IY1 ■*- ™ , ,-. T ri ~^ . ~ 



(2.13) 



2cos(a7r/2) V~°° x ^ 2 ' *~°° dx 2 
For < s < 1, we define 

3-DjM*) + xD~ s u(x) _ -ooI s x u{x) + x I^u(x) 



A_ s/2 u(x) = -- 



2,os(^K) 2cos(f, 



When 1 < a < 2, we have 

To carry out the analysis, we introduce appropriate fractional spaces. 

Definition 2.4. (Left fractional svace\23^ ). We define the following seminorm 

\u\ ja (R) = || ^DSuWvw (2.15) 

and the norm 

i 
H|j£(R) = (Mj£(R) + Nli»(R)) 2 (2-16) 

and let J£(K) denote the closure o/C^°(R) wii/i respect to \\ ■ \\j<*. 

Definition 2.5. (Right fractional space [23] ). We define the following seminorm 

I«Uscr) = II x^S,«IU»w ( 2 - 17 ) 

and i/ie norm 

i 
I|w||j«(k) = (Mj S (r) + IMIl^r)) 2 (2-18) 

and /ei Jfj(R) denote the closure o/C^°(R) with respect to \\ ■ \\j°. 

Definition 2.6. (Symmetric fractional svace l23]/ ) . We define the following 
seminorm 

l«k(R) = \LooD>, x DZu) L2m \" (2.19) 

and t/ie norm 

ll«l|jg(B) = (l«ljf (R) + Nli»(R)) 2 (2-20) 
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and let Jg(R) denote the closure o/Cg°(R) with respect to \\ ■ \\j". 
Using these definitions, we recover the following result, 

Lemma 2.7. (Adjoint vrovertyf3R \2B$) 

(_ 0o />,u) ffi = (u, x J» R (2.21) 

Lemma 2.8. JM/ 

(_ 00 />, x J^u) K = cos(a^)|u|2„ Q(R) = cos(a7r)|u|^ Q(E) (2.22) 

From Lemma I2~7ll2.81 we recover 
Lemma 2.9. For all < s < 1, 

(A_.u,u) a = |u|; U R) = M% (R) (2.23) 



Generally, we consider the problem in a bounded domain instead of R. Hence, we 
restrict the definition to the domain $7 = [a, 6]. 

Definition 2.10. Define the spaces J£ (fl),Jfl (fl),J| (fi)) as £/ie closures of 
Cq°(Q.) under their respective norms. 

For these fractional spaces, we have the following Theorem|18j. 

Theorem 2.11. If -a 2 < -on < 0, then J-^nXorJ^ 1 ^) orJ^ffi)) is 
embedded into J£n 2 (0)(or J^q 2 (&) orJgQ 2 (fl)), and L 2 (f2) is embedded into both of 
them. 

Lemma 2.12. (Fractional Poincare- Fri edrichs) Jjfff For u £ J^ (f2), we ftave 

ll«IU 3 (n) < CMj£ i0 . 

and for u £ J^ (fi), we Ziaue 

IMIi*(n) ^ C\u\j^ Q . 

From the definition of the left and right fractional integral, we obtain the following 
result, 

Lemma 2.13. Suppose the fractional integral is defined in [0,6]. Let g(y) = 
f(b — x), then 

x I£f(x) v = b = x I2g(y) (2.24) 

Lemma 2.14. Suppose u(x) is a smooth function £ fi C R. flh is a discretization 
of the domain with interval width h, Uh{x) is an approximation of u in Pfi. V i, 
Uh(x) £ Ii is a polynomial of degree up to order k, and (u,v)i i — (uh,v)i i ,Vv £ P k . 
k is the degree of the polynomial. Then for — 1 < a $J 0, we have 

||A a/2 u(a:) - A a/2 u h (x)\\ L2 < Ch k+1 
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and for O^n— 1 < a ^ n, k ^ n we have 

|| - (-A)9«(a;) + (-A)%u h (x)\\ L2 < C/i fc+1 -" (2.25) 

where C is a constant independent of h. 

Proof. Case — 1 < a ^ 0, 

First, we consider the approximation error for a fractional integral || a I~ a u{x) — 

First recall that ||u(ir) — Uh(x)\\2 = 0(h k+1 ) from classical approximation theory. 
Suppose x £ Qh,i- We only need to consider r l = a I~ a (0(h k+1 )). 

(x - s)- a ~ l O(h k+1 )ds 



€ 



r(-a) J a 

(b-a)- a O(h k+1 ) 



The case 0<n-l<a<n, recall that |||^ - ^M| < Ch k+1 ~ n and use the 



r(l - a) 

Recalling Lemma f2.13^ the case — 1 < a ^ is proved 

The case < n- I < a < n, rect 
result for case — 1 < a ^ 0, we obtain, 

|| - (-A)*u(a:) + (-A)^u h (x)\\ L2 ^ Ch k+1 - n (2.26) 

This proves the Lemma. D 

Lemma 2.15. (Inverse vroverties U If ) Suppose Vt is a finite element space 
spanned by polynomials up to degree k. For any u>h £ Vh, there exists a positive 
constant C in dependent of Uh and h, such that 

\\d x u h \\ < C/i _1 ||u fe ||, \\u h \\ Th < Ch-^\\u h \\ 

3. LDG scheme for the fractional convection-diffusion equation. Let us 

consider the fractional convection-diffusion equation with 1 < a < 2. From (|2.ip and 
Lemma l2.14[ we know that a direct approximation of the fractional laplacian operator 
is of order k + 1 — n. To obtain a high order discontinuous Galerkin scheme for the 
fractional derivative, we rewrite the fractional derivative as a composite of first order 
derivatives and fractional integral and convert the equation to a low order system. 
Following (|2.14[) . we introduce two variables p, q, and set 

q = A (q _ 2 )/2P 

r d 

p = vs^-u 

ox 

Then, the fractional convection-diffusion problem can be rewritten as 

q = A (q _ 2 )/2P 

Consider the equation in Q = [a, b]. Given a partition a = Xi < xs. < ■ ■ ■ < x K+ i = 6, 
we denote the mesh by h = [x„- i , £„• , i 1, A.xh = .t„- , i — x,- i . 

J •> L J— o J+2 J J J+2 J~ 2 
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We consider the solution in a polynomial space Vh, which is certainly embedded 
into the fractional space according to theorem 12. Ill The piece-wise polynomial space 
Vh on the mesh is defined as, 

V h = {v:v€P k (I j ),x£l j } 

We seek an approximation (uh,Ph, Qh) G Vh to (u,p, q) such that, for any v,w, z £ 
Vh, we have 

' (5s£d if>(aj) ^ + (£/(« fc ), «(*))* = VE(fc,v(*)) /( 

(q h ,w(x)) h = (A {a _ 2)/2 p h ,w(x)) I ^ 3 jx 

(u h (a;,0) ) u(x)) 7i = (u (x), v(x)) u 

To complete the LDG scheme, we introduce some notations and the numerical flux. 
Define 

„.+ i „.- 



u (xj) = lim u(x), Ui 



2 , n-j 

"3 

and the numerical flux as, 

u = h u (u~,u + ), q = h q (q~,q + ), fh = f( u h> u t) 

For the high order derivative part, a good choice is the alternating direction flux [151 
[55] , defined as, 



iij+i = it.,,, j i+ i = 07, i > ° ^i^K -1 



H+ _ 

or 

Ui = u^.i , & = o^i , 1 sc i < if 

For the nonlinear part, /, any monotone flux can be used |28) . 

Applying integration by parts to (|3.1[) . and replacing the fluxes at the interfaces 
by the corresponding numerical fluxes, we obtain, 



((uh)t,v) Ii + yfhv - Veq'hvj \ x '+ 2 - (f(u h ) - y/sq h ,v x ) h = (3.2) 

' l ~2 

{q h ,w(x)) u - (A {a ^ 2 )/2Ph,w(x)) h = (3.3) 

X Z i 
(p h ,z(x))j, ~ \feuhz\ \ 2 +y/e{uh,z x ) I . = (3.4) 

» X i « 

2 

(u/j(x, 0), w(aj)) 7 . - («o(a0) v ( a; ))/ i = ° ( 3 - 5 ) 

Remark 3.1. Originally, the problem is defined in R. However, for numerical 
purposes, we assume there existing a domain tt = [a, b] C R in which u has compact 
support and restrict the problem to this domain Q. As a consequence, we impose 
homogeneous Dirichlet boundary conditions for all u € R\J7. to recover 

. _ A W2 i \ _ qoDx u ( x ) + xD^u(x) _ gD°u(x) + x D%u{x) 
2cos(^J 2cos(^J 
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For the flux at the boundary, we use the flux introduced in |10j . defined as 
U K+ i=u{b,t), q K+ i=q K+ i+jlu K+ i] i 



ui=u(a,t), ?i=?i+^[wi] 

where (3 is a positive constant. 

4. Stability and error estimates. In the following we shall discuss stability 
and accuracy of the proposed scheme, both for the fractional diffusion problem and 
the more general convection-diffusion problem. 

4.1. Stability. In order to carry on analysis of the LDG scheme, we define 



@(u,p,q',v,w,z)= y ^2,{u t ,v) I .dt+ \ ^2(fv~Vsqv)\J+^dt 

.T K _ T K 

~ J2(f{u)-V^q,v x ) I .dt+ y2{q,w(x)) h dt (4.1) 

Jo i=1 Jo l=1 

.T K -T K x - i 

-/ y^{\ a -2)/2P,w{x)) I dt- I ^^feuz\ + 2 dt 

/■T ^ /.T A' t 

+ / 2J (p, z(x)) 7 (it + / y^ -y/e (u, z x ) I . dt — / JC(v,w,z)dt 
Jo i=1 Jo i=1 Jo 

where „Sf contains the boundary term, defined as 

- Zed — 

££{y, w, z) — y/euia, t)zf — u(b, t)v~, i dt + y/su(b, t)z~. t = (4.2) 

2 ft, K + 2 K + 2 

If (it,p, g) is a solution, then .9§(u,p, q; v, w, z) =0 for any (y, w, z). Considering the 
fluxes u ; , i = u~7, i , a,- 1 i = <7 + , i an d the flux at the boundaries we recover, 

! +2 *+2 + 2 *+3 

rT K rT K 

&(u : p,q;v,w,z) = y^(u t ,v) I .dt- / V] (/(u), u a! ) J . di 
Jo i=1 Jo i=1 

+ / ^2(V^q,Vx)j.dt+ y2^/e(u,z x ) I .dt+ ^2(q,w(x)) L 
Jo i=1 Jo i=1 Jo i=1 

,T # /-T K 

-/ y2( A (a-2)/2P,w(x)) I dt+ y^(jp,z{x)) I .dt 

Jo ■_-, * Jo ■_-, 



i=l "" i=l 

/ ^/ i+ |H i+ i*+/ E^ + iHh|* (4-3) 

rT K-l -T 

Jo i=1 Jo 

T /-T r-Q 

V~e(qivi-q K+ ,v KH )dt + ^ —u KH v KH dt 
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Lemma 4.1. Set (v,w,z) — (u,—p,q) in (|4.3|) . and define $(u) = J f(u)du. 
Then the following result holds, 

,%(u,p,q;u,-p,q) = \\u(x,T)\\ 2 - \\u \\ 2 + I {A {a „ 2)/2 p,p)dt + [ v^-( u - fdt 

Jo Jo n 2 

+ j *(u)j - *(«)*+* - (fu)i + {fu) K+h dt 

r T K-l 

Jo j=1 
Proof. Set (v,w,z) = (u,—p,q) in (|4.3|) . and consider the integration by parts 

formula (q,u x )i i + (u,q x )i i = (uq)\ + 2 , to recover at an interface 

x t-i 

2 

A A A-l A'-l 

Y (VeQ,Vx) Zi +'52(Veu,z x ) Ii + Y ^iHi+l + Y ^iWh^ 

i= 1 z— 1 2—1 i— 1 

A x - A'-l K-l 

=^v^(^) i^ + * + Y ^iiW.+i + y ^viM^+i 

i=l i_ 7 i=l i=l 



= — v/e«i «t + Vein-,! 



u. 



J A+I"A+I 

Then 

-T A „T A 



(u,p,q;u,-p,q)= Y {u t ,u) I . dt - / V (f(u),u x ) h dt 

J( > i=l Ja t=l 

rp K t K—\ 

- I Y( A («-v/iP>p)i dt - [ ££+*MI*+i* (4-4) 

J i=l J 2=1 

Define $(u) = J" f(u)du, then 

A A x - A'-l 

2— 1 2—1 X ~^ 2—1 

Combining (|4.4[) and (J4.5I) proves the Lemma. D 

Theorem 4.2. TTie semi-discrete scheme (|3.2[) - (|3.5[) is stable, and \\v,h(x,T)\\ ^ 
||uo(x)|| /or any T > 0. 

Proof. Using the properties of the monotone flux, f(u~,u + ) is a non-decreasing 
function of its first argument, and a non- increasing function of its second argument. 
Hence, we have [^(uh)J - + 1 — fh\uh\ ]+ i. > 0, 1 < j < K — 1. By Galerkin orthogo- 
nality, .^(u h , p h , q h ; u h , ~p h , q h ) = 0, Lemma [O] yields 



jo h [ ~^ 



\\u{x,T)\\ 2 -\\u \\ 2 + (A {a _ 2)/2 p,p)dt+ / v_il(«- ) J dt 



+ 1 *(u)i - *(4+j - (/«)* + (h) K +idt < o 
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Considering the boundary condition and Lemma 12.91 we obtain ||m^(x,T)|| ^ 
|| Wo (x) ||, hence completing the proof. D 

4.2. Error estimates. To estimate the error, we firstly consider fractional dif- 
fusion with the Laplacian operator, i.e., the case with / = 0. For fractional diffusion, 
reduce to 



((uh)t,v)r - {Veqlv)\ + +2 + (\fsq h ,v x ) = (4.6) 

2 

(qh,w(x)) u - (A {a _ 2)/2 p h ,w(x)) h = (4.7) 

_ x i _ 

(ph,z(x))j. - \feu h z\ + +3 +Ve(u h ,z x )j, =0 (4.8) 

s X j » 

2 

(u h (x, 0), w(x)) 7i - («o(x), v(x)) h = (4.9) 

Correspondingly, we have the compact form of the scheme as 

3§(u,p,q;v,w, z) = / V" (u t , v) z . dt + / 'S~\^/e{q,v x ) I . dt + / V] \fe (u, z x ) L dt 
J o i= i Jo l=1 Jo . =1 

+ / y2(Q,w(x))j.dt- y2(A (a _ 2 y 2 p,w(x)) I dt+ y~](jp,z{x)) h dt 
Jo i=1 Jo i=1 Jo f=1 

+ 12 ^ £q t+ iNi+* d *+ / H ^iW 1+ i^ 

/■T" pT rz. o fT 

+ J ^ £q H dt + j — u K +i VK +i dt -J o V-eq KH v KH dt (4.10) 

To prepare for the main result, we first obtain a few central Lemmas. We define 
special projections, ^ 2± and £} into Vt, which satisfy, for each j, 

0> ± u(x)-u{x))v{x)dx = O MveP^ 1 and ^u^x = «(a; ± +i )(4.11) 

(£u(x) - u(x))v(x)dx = Vv e P fc (4.12) 

Denote e« = u — Uh, e p = p — ph, &q = q — Qh, then 

9>- tu = &>- U - u h , &> + e q = &> + q - q h , Be v = £p - p h 
For any (v,w,z) € H^T) x L 2 {Q,T) x L 2 (n,T) , 

3§{u,p,q]v,w,z) — J£{v,w,z) (4-13) 

Hence, <5^(e u , e p , e 9 ; v, w, z) = and we recover 

(3t( ca>- p 6)p c&>+ p ■ cat- p _ 6)p ca>+ p \ 

— (»(' 01>— p p G) P p 02) + p p . CO>~ p 6)p £>}> + p \ 

= S§(^~u -u,3p- p, & + q - q; ,^~e u , Se p , @> + e q ) 
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Substitute {&~u -u,£!p- p, &> + q - q; ^~e u , -£?e p , ^ + e q ) into (|37T0|) . to recover 
the following Lemma, 

Lemma 4.3. For the bilinear form, (|4.10[) . we have 

3S(gp-u -u,£!p- p, ^ + q - q; ,^~e u , -Be p , @> + e q ) 
</ Y,{{^u)t-u u ^e u ) I dt + C T ,a,bh 2k+2 + ^ / y^W^epWJ.dt 

Jo i=1 ' ^%a,b Jo l=1 

where Cx.a.bis independent of h, but may dependent of T and Q. 
Proof. From (|4.10p we have 
^(^~u -u,£p~ p, ^ + q - q; ^~e u , -£e p , 3? + e q ) 



[ S^{{@ > ~u)t-u u & > -e u ) I dt+ I '^^{^ + q-q,{^'e u ) x ) I dt 
Jo i=1 Jo i=1 

r-T K ~T K 

+ / Y, Ve i^ u - u > (^ + e q ) X )j dt - / V {0> + q - q, ^e p ) I dt 

+ V(A (a _ a)/2 (.2j>-p),.2ep) dt+ / y2(£p-p,^ + e g ) r dt 

Jo i=1 ' Jo i=1 



r TK-l _ i-TK-1 

- / 53 Vi(^ + <Z - q)t +i IP-^Ji+i + / E V^(^-u - u)- +k [^ + e,] i+ 1 
J° i=l " J° j=l 



+ | V^ + 9-9)tI^-e+Jrdi + | v|^(^- u _ u )- +4 



IWi df 



7' 



v^(^ + 9 - «)^ + i [^"e-^+idt 

Since (&>-e u ) x £ P*" 1 , (^e^ e P k - l ,(3e p ) x G P k -\j2e p £ P fc , by the proper- 
ties of the projection 3 g± , £?, we recover, 

(<?+q - q, {0>-e u ) x ) u = 0, (&-u - u, {@> + e q ) x ) u = 
(^p-p,^+e q ) Ii =0, (£>p~p,(& + e q ) x ) h =0 
(&~U-u) i+i =0, (& + q-q) t+ i=0. 

We obtain, 

38{^>~u -u,Bp- p, ,9> + q - q; ,^~e u , -£!e p , ,9> + e q ) 



f f2((&-u) t -u u &>-e u ) I .dt- [ ^~e(^ + q-q-) 
Jo i=1 ' Jo 

+ / Y, ( A («-2)/ 2 (^p - p) - (^ + q - g), £e p ) r dt 

J j=l 
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Recalling the projection property and Lemma 12 .141 we obtain \\A/ a _ 2 y 2 (J2p — p) — 
{^ >+ q — q)\\ ^ Ch k+1 . Combining this with Young's inequality and (|2.15j) . we obtain 

3§(^>-u -u,£p- p, 3^ + q - q; @>~e u , -Be p , 3? + e q ) 
</ Y,{{&~u)t-u t ,&-e u ) I dt + C T ,a, b h 2k+2 + 7^— I y)Pcp|| r< 

Jo i=1 ' ^T,a,bJo i=1 

-\{&-e u ) K+ A 2 dt 



2 dt 



/o h 

This proves the lemma. D 



Lemma 4.4. JJQ| / Suppose that for all t > we have 

X 2 (t) + R(t) ^ A(t) + 2 / B(s) X (s)ds, 



JO 

where R, A, B are nonnegative functions. Then, for any T > 0, 

^ X 2 (T) + R(t)^ sup A 1 ' 2 {t) + f B{t)dt 

O^t^T Jo 

Theorem 4.5. Let u be the a sufficiently smooth exact solution to (jl.lj) in VL C M 
mf/i /(u) = 0. Let u^ fee the numerical solution of the semi-discrete LDG scheme 
(|3.2[) - (13.5p . i/ien for small enough h, we have the following error estimates 

\\u-u h \\ ^Ch k+1 

where C is a constant independent of h. 

Proof. From Lemma 14. II and the initial error ||<!3 2 ' _ e u (0)|| = 0, we have 
Vi( ca— p (j) p cai+ p . cat— e ^o) p ca+ p \ 

= \\\^eu(T)\\ 2 + J (A (a _ 2)/2 £e p ,J2e p ) h dt+[ ^\^ e u \\ +l _dt 

Combining this with Lemma 14.31 and (|4.14|) , we have 



1„ rT 



-u(T)\\ 2 + J (A {a _ 2)/2 £e p , £e p ) dt 
</ ((^-u) t -u t ,^-e u )dt + C T ^ b h 2k+2 + -^— f \\£e p \\ 2 dt 

JO ^T,a,b Jo 

Recalling the fractional Poincare-Friedrichs Lemma 12.121 we get 

-||^-e u (r)|| 2 < / ((0>~u) t -u u &>-e u )dt + C T ,a, b h 2k+2 



2 



o 



Using Lemma [4.4I and the error associated with the projection error, proves the the- 
orem. □ 
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For the more general fractional convection-diffusion problem, we introduce a few 
results and then give the error estimate. 

Lemma 4.6. HI] For any piecewise smooth function oj G L 2 (il), on each cell 
boundary point we define 

«(/,«) = «(/,« ,« )-j i |//( _ }| z/ M=Q (4.14) 

where f(u>) = /(u;~,u; + ) is a monotone numerical flux consistent with the given flux 
f. Then k(/,w) is non-negative and bounded for any (lj~,uj + ) G K. Moreover, we 
have 

5l/'(«)l<«(/;w) + a|[«]|, (4.15) 

-i/"(w)[w] <k(/;w) + C»|[w]| 2 . (4.16) 

Define 

if if . if 



if 



.7=1 

Lemma 4.7. fffffl -For J4?(f; u, Uh, v) defined above, we have the following estimate 

K 1 

J2jr j (f;u,u h ,v)^--K(f ] u h )[v} 2 + (C + C4\\v\\ x + h- 1 \\e4l ))\\v\\ 2 

+ (C + C^- 1 ||e u || 2 00 )^ +1 



Combining Theorem 14.51 and Lemma 14.71 we recover the following error estimate 

Theorem 4.8. Let u be the exact solution of (jl.ip . which is sufficiently smooth 
£i!cl and assume f G C 3 . Lei it^ fee i/ie numerical solution of the semi-discrete 
LDG scheme (14.61) - (|4.9|) and denote the corresponding numerical error by e u = u—Uh- 
Vh is the space of piecewise polynomials of degree k Js 1, then for small enough h, we 
have the following error estimates 



\U— Uh\ 



C Ch k+ ^ 



Remark 4.9. Although the order of convergence k + h is obtained it can be im- 
proved if an upwind flux is chosen for f . In this case, optimal order can be recovered. 
We shall demonstrate this in the numerical examples. 

Remark 4.10. For problems with mixed boundary conditions, another scheme 
may be more suitable for implementation. Suppose the problem has a Dirichlet bound- 
ary on the left end and a Neumann boundary at the right end. Let p — v^Tjf > 1 = 



- + £f(u)=A {a _ 2)/2q 
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V^af , then we have 

Si . 

Similarly, we have the following scheme, 

{{uh)t,v) h + hv\u - (f(u h ),v x ) h - (A (a _ 2 y 2 q h ,v(x)) h = (4.17) 

{q h ,w(x)) Li - y/ephMh +Vz{Ph,w x ) u = (4-18) 

{p h ,z(x)) h - ^feu\z\ u + y/e(uh,z x ) u = (4.19) 

(u h {x, 0), v(x)) It - (u (x),v(x)) Ii = (4.20) 

In this scheme, a mixed boundary condition is imposed naturally. However, the anal- 
ysis is more complicated. While computational results indicate excellent behavior and 
optimal convergence, we shall not talk about the theoretical aspect of this scheme. 

5. Numerical examples. In the following, we shall present a few results to 
numerically validate the analysis. 

The discussion so far have focused on the treatment of the spatial dimension with 
a semi-discrete form as 

^ = C h (n h ,t) (5.1) 

where u^ is the vector of unknowns. For the time discretization of the equations, we 
use a fourth order low storage explicit Runge-Kutta(LSERK) method[28] of the form, 



p° = u", 




*e[l,---,5]: 


f k l =a l k l - 1 +AtC h (p z -\t n + c l At) 
\ pW =p( i - 1 )+6 i kW, 


u n+l = p (5) 





where Oj, &j,Cj are coefficients of LSRK method given in [55]. For the explicit methods, 
proper timestep At is needed. In our examples, numerically the condition At ^ 
CAx^ in , (0 < C < 1) is necessary for stability. 

Example 1. As the first example, we consider the fractional diffusion equation 
with the fractional Laplacian, 

^H = -(-Ar/Mx,t)+g(x,t), m [0, 1] x (0,0.5] (52) 

u(x, 0) = Uo(x). 

with the initial condition uo (x) — x e (l— x) 6 , and the source term g(x,t) — e~* (— Uq(x) + (—A)%uo(x)) 
The exact solution is u(x, t) — e" t x e (l — x) e . 

The results are shown in Table |5~T1 and we observe optimal 0(h k+1 ) order of con- 
vergence across 1 < a < 2. 

Example 2. We consider the fractional Burgers' equation 



^ + £ x (^¥ L )=e(-(-Ar/*)u(x,t)+ 9 (x,t), zn [-2, 2] x (O.O.Jjj. 3) 



dt 

u(x, 0) = uq(x) 
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Table 5.1 
Error and order of convergence for solving the fractional Laplacian problem with K elements 
and polynomial order N. 





a = 1.1 


K 


10 


20 




30 




40 




N 


IM|a 


IM|a 


Order 


IMIa 


Order 


IMIa 


Order 




K= 10 


K= 20 


order 


K= 30 


order 


K= 40 


order 


1 


1.64e-06 


3.98e-07 


2.04 


1.73e-07 


2.05 


9.60e-08 


2.05 


2 


1.24e-07 


1.69e-08 


2.87 


5.13e-09 


2.94 


2.18e-09 


2.97 


3 


1.14e-08 


7.23e-10 


3.98 


1.41e-10 


4.03 


4.44e-ll 


4.03 




a = 1.3 


K 


10 


20 




30 




40 




N 


IM|a 


||e„|| 2 


Order 


IMIa 


Order 


IMIa 


Order 




K= 10 


K= 20 


order 


K= 30 


order 


K= 40 


order 


1 


1.45e-06 


3.46e-07 


2.07 


1.50e-07 


2.06 


8.33e-08 


2.05 


2 


1.19e-07 


1.55e-08 


2.94 


4.64e-09 


2.98 


1.96e-09 


2.99 


3 


1.04e-08 


6.48e-10 


4.00 


1.26e-10 


4.03 


3.97e-ll 


4.02 




a= 1.5 


K 


10 


20 




30 




40 




N 


IM||a 


IM|a 


Order 


IMIa 


Order 


1Mb 


Order 


1 


1.35e-06 


3.24e-07 


2.06 


1.42e-07 


2.04 


7.90e-08 


2.03 


2 


1.22e-07 


1.51e-08 


3.01 


4.51e-09 


2.99 


1.90e-09 


2.99 


3 


1.04e-08 


6.28e-10 


4.05 


1.22e-10 


4.04 


3.85e-ll 


4.02 




a= 1.8 


K 


10 


20 




30 




40 




N 


IMIa 


1Mb 


Order 


IMIa 


Order 


1Mb 


Order 


1 


1.28e-06 


3.11e-07 


2.04 


1.38e-07 


2.01 


7.74e-08 


2.01 


2 


1.40e-07 


1.51e-08 


3.21 


4.48e-09 


3.00 


1.89e-09 


3.00 


3 


1.44e-08 


6.72e-10 


4.42 


1.23e-10 


4.19 


3.83e-ll 


4.05 



with initial the condition 



u (x) = 



(l-x 2 ) 4 /10 




-1 <a;< 1 

otherwise 



The source term g(x,t) = e * (— uq(x) + e t Uo(x)u' (x) + e(— A) 2 uq(x)) ,e = 1 is 
derived to yield an exact solution as 



u(x,t) 



e-*(l-a; 2 ) 4 /10 




-1 < x < 1 

otherwise 



The results are shown in Table [5T^1 and we observe optimal 0(h k+1 ) order of conver- 
gence across 1 < a < 2. 



Example 3 Let us finally consider the fractional Burgers' equation with a dis- 
continuous initial condition, 

u (x) = { 3 - 1 ^ x ^° 
I otherwise 

We consider the case in (|5.3|) with e — 0.04, g(x,t) = 0. The increasing dissipative 
effect of increasing a is clear from the sequence of figures, Fig. (|5.1I) - (|5.4[) . with the 
latter case of a — 2 being the classic case. 
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Table 5.2 
Error and order of convergence for solving the fractional Burgers' equation with K elements 
and polynomial order N. 





a= 1.01 


K 


10 


20 




30 




40 




N 


IM|a 


1Mb 


Order 


1Mb 


Order 


1Mb 


Order 




K= 10 


K= 20 


order 


K= 30 


order 


K= 40 


order 


1 


1.10e-03 


2.81e-04 


1.97 


1.24e-04 


2.03 


6.90e-05 


2.03 


2 


6.53e-05 


1.00e-05 


2.70 


3.09e-06 


2.90 


1.33e-06 


2.94 


3 


5.94e-06 


4.00e-07 


3.89 


8.05e-08 


3.96 


2.58e-08 


3.95 




a = 1.5 


K 


10 


20 




30 




40 




N 


IMIa 


1Mb 


Order 


IMIa 


Order 


IMIa 


Order 


1 


8.89e-04 


2.15e-04 


2.05 


9.45e-05 


2.03 


5.28e-05 


2.02 


2 


6.71e-05 


8.62e-06 


2.96 


2.57e-06 


2.99 


1.09e-06 


2.99 


3 


4.91e-06 


3.34e-07 


3.88 


6.80e-08 


3.93 


2.16e-08 


3.99 




a= 1.8 


K 


10 


20 




30 




40 




N 


IMIa 


1Mb 


Order 


IMIa 


Order 


IMIa 


Order 


1 


8.43e-04 


2.09e-04 


2.01 


9.25e-05 


2.01 


5.20e-05 


2.00 


2 


6.78e-05 


8.59e-06 


2.98 


2.56e-06 


2.99 


1.08e-06 


2.99 


3 


4.80e-06 


3.32e-07 


3.85 


6.84e-08 


3.90 


2.20e-08 


3.94 





Fig. 5.1. Solution of the fractional Fig. 5.2. Solution of the fractional 

Burgers's equation with e = 0.04, a = 1.005 Burgers's equation with e = 0.04, a = 1.1 



6. Concluding remarks. We propose a discontinuous Galerkin method for 
convection-diffusion problems in which the a fractional diffusion is expressed through 
a fractional Laplacian. To obtain a high order accuracy, we rewrite the fractional 
Lapacian as a composite of first order derivatives and integrals and transform the 
problem to a low order system. We consider the equation in a domain $7 with ho- 
mogeneous boundary conditions. A local discontinuous Galerkin method is proposed 
and stability and error estimations are derived, establishing optimal convergence. In 
the numerical examples, the analysis is confirmed for different values of a. 
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Fig. 5.3. Solution of the fractional Fig. 5.4. Solution of the fractional 

Burgers's equation with e = 0.04, a = 1.5 Burgers's equation with e = 0.04, a = 2.0 
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